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SUMMARY 

The  Finite  Element  and  Boundary  Element  Methods  are 
described  with  their  essential  features  illustrated  using  an  example 
of  a  boundary  value  problem  for  a  harmonic  function.  Analysis  of 
the  methodical  errors  is  then  carried  out.  This  is  followed  by  a 
consideration  of  the  relative  computational  advantages  of  the  two 
methods. 
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1 .  INTRODUCTION 


In  conducting  studies  of  aircraft  flight  dynamic 
behaviour  it  is  often  necessary  to  estimate  the  aerodynamic 
characteristics.  Methods  such  as  those  embodied  in  the  USAF 
Stability  and  Control  DATCOM  and  in  the  ESDU  data  sheets  are 
available,  but  they  are  rather  empirical  with  all  their  inherent 
limitations.  Recent  developments  in  computational  fluid  dynamics 
offer  the  prospect  of  estimating  aerodynamic  characteristics 
directly  from  physical  principles  at  a  cost  acceptable  for 
engineering  purposes. 

As  a  first  step  towards  adopting  such  methods,  an 
examination  of  the  relative  merits  of  the  Finite  Element  and  the 
Boundary  Element  methods  is  made.  The  purpose  of  this  Memorandum 
is  to  provide  a  concise  and  intelligible  summary  of  the  two  methods, 
with  sufficient  emphasis  on  their  essential  features.  Such  a 
summary  is  thought  to  be  desirable  in  the  face  of  a  proliferation 
of  bewildering  and  often  fragmented  literature  on  the  subject  (for 
example,  see  the  list  of  references  of  [1]  to  [4]).  The  example 
chosen  to  illustrate  the  methods  is  the  boundary  value  problem  for 
a  harmonic  function.  This  problem  is  considered  typical  as  such 
familiar  problems  as  elastostatics  and  potential  flows  can  all  be 
reduced  to  the  problem  of  determining  harmonic  functions. 

The  analysis  of  methodical  errors  is  then  presented. 

It  is  surprising  that  such  essential  analysis  is  not  readily 
accessible  in  spite  of  the  large  amount  of  literature  on  the  two 
methods.  The  claim  that  the  Boundary  Element  Method  is  computation¬ 
ally  advantageous  compared  with  the  Finite  Element  Method  is  also 
examined.  This  claim  is  hard  to  substantiate  when  the  test  problem 
concerns  a  three  dimensional  region. 

Finally,  an  appendix  proves  the  convergence  of  the 
Gauss-Seidel  iteration  method  for  positive  definite  Hermitian  matrices. 
Such  matrices  often  arise (s)  from  the  use  of  the  Finite  Element  Method 
involving  complex-valued  functions. 


2.  FORMULATION  OF  THE  FINITE  ELEMENT  METHOD 

This  method  depends  of  the  existence  of  a  functional 
I,  the  extremal  value  of  which  gives  the  required  solution  to  our 
particular  problem  [2,3,4].  To  fix  ideas,  we  consider  the  problem 
of  finding  a  real-valued  function  $(x)  which  has  continuous  second 
order  derivatives  in  a  region  f 2,  <f>eC2(ft),  such  that 


«  0  in  the  region  ft 


..(1) 


subjected  to  the  boundary  conditions 


4>  =  u  on  , 


..(2) 


[2] 


=  p 

In 


on  r. 


..(3) 


where  the  two  disjoint  surfaces  r,  and  form  the  boundary  T  of  the 
region  ft.  (The  region  ft  and  its  Boundary  T  are  such  that  all  class¬ 
ical  theorems  of  calculus  are  applicable).  It  is  straightforward  to 
use  the  divergence  theorem  to  prove  that  the  above  problem  has  only 
one  solution  0(x).  When  F  =  the  above  boundary  value  problem 
becomes  the  Dirichlet  problem  and  when  T  =  Fg,  the  Neumann  problem 
(only  in  this  latter  case  0(x)  is  determined  only  up  to  a  constant). 
Here  we  take  for  granted  that  there  exists  a  solution  0(x)  to  the 
above  boundary  value  problem. 

Let  us  form  the  functional 


i(< i>)  =  N|2  -  2 


..(A) 


where  p  -  p^  on  F ^  and  p  -  0  on  r.  .  Minimization  of  this  functional 
I  with  respect  to  f(x)  such  that  |  =  n  on  gives 


60  Vz0 


for  any  arbitrary  incremental  function  60  £  C2  (Q)  satisfying  50  *  0 
on  T  .  Hence,  instead  of  solving  for  0  satisfying  the  system  of 
equations  ((1),  (2),  (3)}  we  can  minimize  (4)  subjected  to  the 
condition  (2).  This  is  the  basis  for  the  Finite  Element  Method. 

Thus  the  Finite  Element  Method  replaces  the  function 
0  e  Cz(ft)  by  an  approximation  function  <I>(x)  which  is  determined  by 
a  finite  number  of  its  nodal  values  0(x^),  m  -  1,  N  inside  ft.  These 
nodal  values  are  determined  by  the  minimization  of  J(0).  The  non- 
nodal  values  0(x)  are  determined  from  the  N  nodal  values  so  obtained. 
In  all  the  following,  the  function  0(x)  in  each  elemental  volume  will 
be  a  linear  interpolation  of  its  nodal  values  at  the  vertices  of  the 
element,  as  a  result,  the  function  0(x)  satisfies  V2O  =  0  within  each 
elemental  volume  but  not  across  inter-element  boundaries  (the  higher 
order  interpolation  of  4(x)  does  not  automatically  satisfy  the 
condition  Vz0  «  6  within  each  element  and  additional  equations  are 
imposed  as  a  result) . 

The  functional  J(0)  is  now  discretized  into 


ft-  U  V. 


r-(Yi  Uy2) 


2 

2  E 
i=l 


..(6) 


as  in  Figure  1. 


[3] 


Since  $(x)  within  each  elemental  volume  can  not  be 
affected  by  any  nodal  values  outside  the  element,  we  have  the 
equations  determining  as  in  the  following: 


V 


For  a  node  inside  fi,  i.e.  not  on  the  boundary 


7 

0  =  Z 

._j  V  geometry  of  (a 


f real  constant 

«  V 


depending  only  on  the  A  ^ 
liCL2’  a3J  a4J  aSJ  a6*  a7^  )  ai 


..(7) 


For  a  node  on  the  boundary 


v  /real  constant  depending  only  on  A 

P“3  =  i-1,2,3,4  \the  8eometry  of  ) 


a . 

i 


..(8) 


V 


For  typical  nodes  a^,  belonging  to  the  boundary 


..(9) 


All  equations  of  the  kind  (7),  (8),  (9)  form  a  large 
system  of  n  equations  in  n  unknowns.  The  function  4>(?)  Is  considered 
known  when  its  approximation  4>(a:)  can  be  shown  to  be  reasonably  close 
to  it. 


The  right  hand  sides  of  equations  (7)  and  (8)  result 

from  the  differentiation  of  a  homogeneous  second  degree  polynomial 

in  <J  .  Therefore,  we  have 

a. 


N 

Z 

i-1 


o.  3$ 

%  a . 


|v<t>| 2  >  o  . 


Hence  the  matrix  A  of  coefficients  in  the  equations  (7),  (8)  and  (9) 
satisfies 


a  -o . 

t  3 


..(10) 


and  is  positive  definite  (it  is  easy  to  show  that  A  is  also  symmetric). 
The  Gauss-Seidel  iteration  method  can  thus  be  used  to  advantage  in 
solving  for  the  unknowns  of  the  system  {  (7) , (8) , (9)} ,  see  the 
appendix  at  the  end  of  the  memo. 


[4] 


The  calculations  in  this  section  also  hold  if  the 
function  $(£)  is  complex-valued.  In  this  case  the  last  integral 
of  equation  (4)  changes  to  Real  (f  <p  p)  where  p  is  the  complex 
conjugate  of  p,  the  matrix  A  of  equation  (10)  becomes  a  positive 
definite  Hermitian  matrix,  also  equations  (5)  and  (6)  change 
slightly  but  this  does  not  alter  our  final  system  { (7) , (8) , (9) }  . 


3.  FORMULATION  OF  THE  BOUNDARY  ELEMENT  METHOD 


To  solve  for  the  solution  of  the  same  system  {(1),(2), 
(3)}  the  Boundary  Element  Method  starts  with  the  identity 


(V2$)  y  =  0  ..(11) 

ft 


where  y  is  any  C2(ft)  function.  Using  Green's  theorem  to  integrate 
the  above  integral  by  parts,  we  have 


0  (V2z/) 
ft 


9m 


9  <;> 

9m 


=  0 


Jr 


..(12) 


Letting  V2y  =  -6($  - %*)  where  x*  is  on  the  boundary  T 
of  the  region  ft,  we  have  the  following  singular  equation 


Y  4>(x*)  + 


¥  9m 


r  -  six*) 


y&.O 
'  9m 


T -S(x*) 


..(13) 


where  Six*)  is  a  very  small  sphere  centered  on  x* e T  and  y  is  a 
constant  equal  to  the  fraction  of  the  surface  of  Six*)  contained 
in  ft;  when  T  is  smooth  this  ratio  Y  is  equal  to  one  half. 

The  function  4>eC2(ft)  is  now  approximated  by  a 
function  'F(x)  which  is  determined  by  its  nodal  values  Vix  ), 
m  *  I,  M  on  the  boundary  T.  The  function  ¥  is  harmonic  and  its 
value  on  each  elemental  boundary  element  can  be  taken  to  be 
constant  or  the  interpolation  between  the  nodal  values  at  the 
vertices  of  the  (boundary)  element.  The  equations  determining 
the  nodal  values  of  V  are  thus 


Y  ¥(x.)  +  E 

**  **  J 


3i /*• 
9m 


M 


ri 


S(x{) 


-  E 

0=1  J 


9T  . 
Vi! n  =  0 


r  .-six.) 
o 


i»  •  •  ■ 


..(14) 


where  y.ix)  is  given  by 
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hi  <*) 


J_  1 

4R  j  x  -  x . I 

1  -  ~t 


i 


1: 


.M. 


..(15) 


In  any  case  the  equations  (14)  have  only  M  unknowns 
since  either  ’{'(x.)  or  j3*f  v  is  known  for  a  given  node. 

3 n 

For  nodes  x.  on  I*,,  we  have 

-V  t  1 


G.  .  (x.)  =  H.  .  V(x.)  , 

ij  3 n  -t  ic  ’ 


1  1  .,  a  .  a  A/  ^  a 


For  nodes  x.  on  f. 


H.  .  ¥(x.)  -  J. .  #  (x.) 
ij  tj  3n 


i  *  A^  +  lt . . .  jM. 


••(16) 


The  unknows  on  the  right  hand  side  of  equations  (16) 
are  then  solved  and  the  value  of  H'(r)  at  any  interior  point  a  of 
is  easily  worked  out  using  equation  (12)  with  V2y  =  ~S(x  -  a) . 

Note  that  the  approximation  function  'l'(x)  so  derived 
satisfies  equation  (1)  everywhere  but  satisfies  the  boundary  conditions 
(2)  and  (3)  only  at  A.'  points.  In  other  words,  equation  (11)  is  always 
satisfied  inside  P  with  H  (x)  while  equation  (13)  on  the  boundary  T, 
which  should  also  be  satisfied  identically  for  every  singular  function 
z  (x)  =  l/|x-x*|,  x*eT,  is  satisfied  for  only  A.'  functions  y^(x) 
defined  by  equation  (15).  For  this  reason  the  method  has  been 
considered  by  some  authors  a  weighted  residual  method,  meaning  that 
equation  (13)  is  satisfied  only  for  a  finite  number  ( M ,  in  this  case) 
of  weighting  functions  y^(x). 

It  is  now  obvious  that  the  Boundary  Element  Method 
depends  on  the  use  of  Green's  identity 

•  ► 

«.V=i>  -  rVM  ..(U, 

although  the  usual  approach  using  the  weighted  residual  method  often 
disguises  this  under  the  process  of  partial  integrations.  Another 
interesting  feature  of  the  Boundary  Element  Method  is  the  applicability 
of  the  Maximal  Value  Theorem  which  states  that  a  harmonic  function 
attains  its  maximal  values  only  on  the  boundary  [5].  The  value  of 
the  methodical  error  |$(x)  -  H  (x)  |  at  an  interior  point  of  (2  is  thus 
always  less  than  such  value  at  some  point  on  the  boundary  T  and  the 
latter  can  be  reduced  with  a  finer  discretization  of  T. 


[6] 


4.  ERROR  ANALYSIS  FOR  THE  BOUNDARY  AND  FINITE  ELEMENT  METHODS 

The  determination  of  the  maximum  error  | (a:)  -  Hr)| 
or  |4>(x)  -  $(r) |  in  Boundary  and  Finite  Element  Method  is  not  at  all 
simple.  Although  the  Maximal  Value  Theorem  gives  an  upper  limit  for 
the  error  function  in  the  Boundary  Element  Method  we  are  unlikely  to 
know  the  limit  unless  ■  T. 

4.1  Boundary  Element  Method 


Let  us  first  examine  the  simpler  case  of  the  Boundary 
Element  Method.  We  can  use  Green's  identity  to  write 


4JI  [<j>(x)  -  Ux)] 


<4>  -  V)  (^y)  +  l  (U  -  fQ 

n  T~\ 


..(18) 


we  can  say  that  the  point  error  at  each  point  depends  mostly  on  the 
error  on  the  boundary  points  closest  to  it  and  the  variation  is  more 
prominent  with  (<  -  ¥)  than  with  ( 3d>  8H'\. 


rith  /jhj>  8H'  \ . 

~  J 


To  remove  all  the  unknown  quantities  such  as  (tp  -  T) 
on  F  and  / \  on  T  we  have  to  use  a  Green's  function  G(y)  satis¬ 
fying  \dr.  )  1 


V2  G(y)  =  o(y  -  x)  , 


;<x)  =  0  on  L 


and  igp  =  0  on 


..(19) 


A  new  relation  is  then  arrived  at 


4Ti  [*<*)  -  Hr)] 


(4>-W) 


9G  ,  r  (*$L 
dZ+  G\dn  dn) 


..(20) 


Obviously,  a  new  function  G(x )  has  to  be  calculated  but  the  error 
(<{)(x)  -  H'(x)]  can  be  determined  without  any  ambiguity.  (The  pole 
y  of  the  function  G  has  been  assumed  to  be  in  the  interior  of  1). 

If  T  is  smooth  and  j/tF  then  at  least  one  of  the  integrals  on  the 
right  hand  side  is  Improper  and  the  left  hand  side  must  correspond¬ 
ingly  be  divided  by  two) .  Another  inequality  will  also  be  given  in 
the  next  subsection. 


[7] 


4.2  Finite  Element  Method 


The  estimation  of  error  for  the  Finite  Element  Method 
is  more  complicated  as  we  have  here  in  addition  to  the  interpolation 
error  on  the  boundary  another  kind  of  error  arising  from  the  dis¬ 
cretization  of  the  function  C(x).  This  makes  its  gradient  VC  vary 
in  a  stepwise  manner  within  ft  and  the  jumps  occur  at  inter-element 
boundaries .  The  Maximal  Value  Theorem  is  not  applicable  in  this 
case  as  $  is  not  harmonic  over  the  region  ft. 

A  natural  way  to  measure  the  closeness  of  C(x)  to 
C(x)  is  to  use  the  difference 


1(C)  -  Kc) 


(  |VCj2  -  |V$|2)  -  2 


(C  -  4>)  p 


Hence 


1(C)  -  1(c)  =  2 


VC  .  V(C  -  <j>)  - 


[  V(C  -  C)  ]: 


..(21) 


<■> 


a  -Op. 


..(22) 


Hence  to  have  more  accurate  results  we  must  make 
V(C  -  C)  small  especially  when  VC  is  large  and  similarly  for  (C  -  <J>) 
on  the  boundary,  it  must  be  small  where  p  is  large. 

Similarly  to  the  case  of  Boundary  Element  Method,  the 
error  in  Finite  Element  Method  can  be  written  down  as 


411  (C(x)  -  C  (r)  ] 


••(23) 


where  G  is  the  Green's  function  defined  by  equation  (19)  and  S  is  a 
surface  comprised  of  the  double  layers  enclosing  all  inter-element 
boundaries  with  the  normals  pointing  away  from  the  elemental  volumes. 
The  equation  holds  for  any  point  x  interior  to  ft.  If  *  is  right  on 
the  boundary  T  then  at  least  one  of  the  first  two  integrals  on  the 
right  hand  side  of  (23)  is  improper  and  also  the  left  hand  side  must 
be  divided  by  two.  The  last  integral  in  equation  (23)  is  the  error 
caused  by  the  stepwise  variation  of  VC.  Therefore,  the  errors 
(C  -  C)  and  ( ci  _  H  N  on  Tj  and  respectively,  as  well  as  the 
\en  Sn  y 

jump  in  VC  across  inter-element  boundaries  must  all  be  reduced  if  we 
want  to  reduce  the  magnitude  of  the  left  hand  side  of  (23). 


[8] 


An  often  suggested  measurement  for  the  global  performance 
(i.e.  performance  over  the  whole  region  ft  in  an  average  manner)  of  a 
Finite  Element  scheme  is  the  functional  J  defined  by 


-  *) 


|  *-*  |2  + 


9£  _ 

in  in 


In  the  following  we  will  examine  the  usefulness  of  this  functional. 


From  equations  (23)  and  (24)  we  have  an  inequality 


|<fr(x)  -  i(x) 


fx'H  nr  *  |  w1]1 

'•Jr.  Jr„  ■> 


which  gives  the  upper  bound  for  the  point  error  [<f>fcc)  -  $0*:)]  in 
terms  of  J  and  the  Green's  function  defined  by  eo  on  (19).  It 
is  obvious  that  even  if  we  are  only  concerned  wi  the  resulting 

values  of  £(.r)  and  C(r)  on  the  boundary  T  for  a  ticular  Finite 

Element  scheme,  the  error  originating  from  inte  lenient  discontin¬ 
uities  in  Vi  still  has  to  be  accounted  for.  He  .e  .e  performance 
of  the  scheme  can  not  be  adequately  measured  by  ,  use  of  J  alone. 


On  the  other  hand,  if  we  omit  the  first  term  on  the 
right  hand  side  of  equation  (25)  and  replace  <t>(x)  by  ¥(x)  we  have 
an  inequality  for  the  Boundary  Element  Method.  'Hence  J~is  a  good 
measurement  of  the  global  performance  of  any  Boundary  Element  scheme. 

«/  can  also  be  related  to  the  quantity  [/($)  -  J(< J>)] 
by  derivation  from  equations  (21)  and  (24) .  We  get 


iw  -  KO  l  < 


($  -  <fr) 


Is 


.(26) 


The  first  integral  is  the  effect  of  the  jump  in  Vi 
across  inter-element  boundaries  and  the  last  integral  becomes  negligible 
when  $  is  close  to  i. 


In  a  similar  way  we  also  have 


|  I(i)  -  1(0  |  < 


j  I  *  -  <t>  I2  + 

Js  Jr 


i  *  -  ♦  I2  (  * 


+  !£  +  !£  + 
3n  3n  Tin 

5  Jr. 


9$  n  *  t 

3n  "  p  {  •  *  * 1 

r„  J 


which  shows  the  effect  of  |<t  -  <j>|  on  the  difference  J(<2>)  -  I(4>) . 


5.  RELATIVE  COMPUTATIONAL  ADVANTAGES  OF  FINITE  AND  BOUNDARY 
ELEMENT  METHOD  AS  FOR  THE  TEST  PROBLEM 


The  case  where  the  region  ft  is  one  dimensional  is 
much  too  simple  to  have  any  significance.  Therefore,  only  the 
cases  where  ft  is  two  or  three-dimensional  are  considered. 

First  when  ft  is  two-dimensional  the  FEM  (Finite 
Element  Method)  has  n2  nodes  and  requires  n 2  memory  locations.  In 
contrast,  the  BEM  (Boundary  Element  Method)  requires  only  n  nodes 
but  also  n2  memory  locations.  If  the  Gauss-Seidel  iterative  method 
is  used  for  the  banded  matrix  of  the  FEM  and  the  standard  Doolittle 
or  Gaussian  method  is  used  for  the  full  natrix  of  the  BEM  then  their 
computing  time  will  be  proportional  to  n*  and  n 3  respectively.  Thus 
the  BEM  is  advantageous  in  this  case. 

On  the  other  hand,  when  ft  is  three-dimensional,  the 
FEM  and  BEM  requre(s)  n3  and  n 2  nodes  respectively.  Memory  locations 
required  are  of  order  n3  and  nu  respectively  and  the  computing  times 
are  the  same.  Thus,  the  FEM  is  advantageous  in  this  case. 


6.  FINAL  REMARKS 


Besides  those  points  considered  in  the  previous  section, 
others  which  favours  one  method  over  the  other  are: 

(i)  The  gradient  of  <p(x)  is  approximated  by  step-like  function 
VS  in  FEM  compared  to  a  harmonic  function  in  BEM. 

(ii)  The  setting  up  of  a  BEM  program  requires  some  elaborate 
calculation  for  each  singular  component  of  the  surface 
integral.  This  is  in  direct  contrast  with  the  straight¬ 
forward  setting  up  procedure  of  FEM.  Moreover,  the  values 
of  the  integrals  in  BEM  are  fairly  sensitive  to  the 
accuracy  of  the  integration  near  to  their  singularities. 

This  creates  yet  another  kind  of  error  which  has  to  be 
considered. 

(iii)  For  problems  with  a  non-homogenous  body  the  BEM  requires 
a  subdivision  of  the  body  while  FEM  has  this  feature 
inherent  in  it. 

(iv)  With  any  subdivision  of  a  body  such  as  in  (iii)  or  as  for 
a  long  a  narrow  region  (for  example,  see  [I],  p.  183)  the 
BEM  becomes  a  hybrid  between  itself  and  the  FEM.  Mixed 
characters  can  thus  be  expected  in  this  case  :  For  example, 
the  coefficient  matrix  becomes  slightly  banded  and  the 
function  'f(x)  is  not  harmonic  at  all  interior  point  of  the 
region  ft. 


The  truncation  error  in  the  Gauss-Seidel  iteration  for 
the  FEM  is  usually  of  greater  order  of  magnitude  than 
the  round  off  error.  The  latter  also  does  not  improve 
much  with  the  use  of  internal  double  precision  in  its 
equation  solving  subroutine. 
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APPENDIX 


Convergence  of  Gauss-Seidel  Method  for  Positive  Definite 
Hennitian  Matrices 


To  make  the  memo  self-contained,  a  proof  of  convergence 
for  the  Gauss-Seidel  method  as  applied  to  a  positive  definite 
Hermitian  matrix  is  given  here.  This  proof  is  a  generalization 
of  the  one  given  in  [7],  p.  445  for  a  real  symmetric  matrix.  The 
convergence  is  independent  of  the  initial  vector. 

Consider  the  system  of  equation 


A  x  =  b 


where  A  is  a  positive  definite  Hermitian  matrix. 
Decompose  A  into 


A  =  L  +  D  +  U 


where  L  and  V  are  strictly  lower  and  upper  complex  matrice  and  D 
a  diagonal  positive  definite  real  matrix.  The  Gauss-Seidel  iteration 
method  then  gives 


x ...  =  -(Z7  +  L)"2  U  x.  +  ( D  +  L )‘J  b 


Hence , 


The  sequence  {x.}  converges  if  all  the  eigenvalues 
of  the  matrix  B  m  -(D  +  f-)”1  i/'feave  their  magnitudes  smaller  than 
unity. 

Since  A  is  Hermitian  we  have 

-I  —T 

B  -  -  (D  +  L)  L 

where  the  overbar  denotes  the  complex  conjugate.  Let  A  and  v  be  an 
eigenvalue  and  eigenvector  of  the  matrix  B.  Ve  will  prove  that 

IM  < 


A. 2 


Since, 


tT 


L  v  ■  \{D  4  L)v  , 


or 


( 2  4  A)  t>  ( D  4  Z,)y  ■  v  A  v  >  0  , 

A  is  different  from  -2. 

Taking  the  conjugate  transpose  of  the  above  equation 

we  have 

(2  4  A)  V  (D  4  L)  v  -  (2  4  A)  y  (£>  4  LT)  y 

-  (2  4  A)  y  0  tf  4  A(2  4  A)  V  {D  4  L)  V, 


or 


(2  -  |A|a)  v  (0  4  I)  r  «  (2  4  A)  y  Z)  y 


So, 


(2  -  |  A | 2)  3  A  y  -  (2  4  A)  (2  4  l)  v  D  v  . 


The  right  hand  side  is  greater  than  zero  since  D  is 
positive  definite.  Thus  |A|  <2,  which  proves  the  convergence'of 
the  iteration  scheme. 
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